function [an, bn] = mie_theory_fn(m, x, nmax)

% INPUT: 
% complex refractive index m=m'+im"
% Size parameter x = k0*a, where k0 = wave number in the ambient medium, a = sphere radius
% Number of terms in expansion nmax

% OUTPUT:
% Mie coefficients an, bn

% Adapted from a code by C. M�tzler, June 2002
% http://arrc.ou.edu/~rockee/NRA_2007_website/Mie-scattering-Matlab.pdf

% create n array
n   = (1:nmax);  
nu  = n + 0.5;
y   = m.*x; 
m2  = m.*m;

sqx = sqrt(0.5*pi./x); 
sqy = sqrt(0.5*pi./y);

% bessel function quantities
bx  = besselj(nu, x).*sqx;
by  = besselj(nu, y).*sqy;
yx  = bessely(nu, x).*sqx;

hx  = bx + 1i*yx;
b1x = [sin(x)/x,  bx(1:nmax-1)];
b1y = [sin(y)/y,  by(1:nmax-1)];
y1x = [-cos(x)/x, yx(1:nmax-1)];
h1x = b1x    + 1i*y1x;
ax  = x.*b1x - n.*bx;
ay  = y.*b1y - n.*by;
ahx = x.*h1x - n.*hx;

% finally the coefficients themselves
an  = (m2.*by.*ax - bx.*ay)./(m2.*by.*ahx - hx.*ay);
bn  = (by.*ax     - bx.*ay)./(by.*ahx     - hx.*ay);

return

